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1.  Introduction 


Additive  manufacturing  (AM)  is  a  production  method  that  involves  gradual,  layer- 
by-layer  building  of  material,  rather  than  traditional  methods  such  as  milling  ma¬ 
terial  from  a  blank  (subtractive  manufacturing)  or  casting  molten  material  into  a 
mold.  There  are  several  AM  methods  for  a  wide  range  of  materials,  such  as  ex¬ 
trusion  deposition,  light  polymerization,  selective  laser  melting,  and  selective  laser 
sintering.  Each  of  these  methods  can  be  used  with  various  materials  such  as  poly¬ 
mers,  ceramics,  and  metal  alloys. 

While  AM  has  greatly  expanded  the  design  space — allowing  the  production  of  pre¬ 
viously  unmanufacturable  topologically  optimized  structures — constraints  remain. 
One  constraint,  for  example,  is  simply  the  minimum  allowable  feature  size  of  a 
manufacturing  technique.  In  other  words,  each  technique  has  a  lower  limit  on  the 
size  of  fine  features  that  may  be  manufactured.  One  method  of  imposing  a  minimum 
feature  size  is  to  use  a  Heaviside  projection  method.1  While  the  minimum  feature 
size  is  a  necessary  constraint  for  any  AM  process  (a  feature  size  smaller  than  the 
minimum  printable  feature  size  is  simply  impossible),  there  are  some  constraints 
that  impact  the  success  rate  of  AM  builds.  One  such  AM  constraint  specifically  re¬ 
lated  to  powder  bed  methods  is  part  failures  due  to  overhanging  material  without  the 
inclusion  of  sacrificial  supports.2  Powder  bed  method  part  failures  typically  occur 
because  of  unsupported,  overhanging  features  that  curl  or  warp  under  thermal  load 
and  are  subsequently  struck  by  the  recoater  blade/roller.  Support  structures  act  to 
wick  heat  away  and  provide  structural  support  to  prevent  thermally  induced  curling. 

One  approach  to  correcting  the  overhang  issue  is  to  incorporate  strict  overhang  con¬ 
straints  on  the  optimized  geometry.3  Alternatively,  a  physics-based  method  may  be 
used  wherein  a  process  model  is  incorporated  into  the  optimization  method  and 
some  parameter  of  that  process  is  optimized.  Here,  we  attempt  to  minimize  thermal 
distortions  induced  by  the  manufacturing  process,  which  should  reduce  the  need  for 
sacrificial  supports.  Because  a  full-resolution  process  model  would  be  too  compu¬ 
tationally  expensive  to  use  in  an  optimization  loop,  we  propose  2  surrogate  models 
for  the  manufacturing  process.  The  first  is  a  basic  thermomechanical  model  that 
assumes  that  the  built  part  is  heated  to  a  uniform  initial  temperature  (T0)  and  cools 
to  a  final  temperature  (Tf).  The  second  is  a  thermomechanical  element-birth  model 
in  which  elements  are  activated  sequentially,  and  activated  elements  have  a  given 
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initial  temperature.  The  activated  elements  then  cool  according  to  a  given  thermal 
history,  eventually  reaching  a  final  temperature  as  before.  This  surrogate  model 
simulates  the  AM  process  more  closely  because  in  direct  metal  laser  sintering  (and 
related  techniques),  features  are  birthed  through  the  laser  heat  source  hitting  the 
powder,  selectively  melting  a  small  section  of  the  powder  bed.  This  molten  material 
then  cools  to  solidify  as  a  solid  metallic  feature  within  the  part  topology. 

2.  Problem  Definition 

In  this  section,  overviews  of  the  compliance  minimization  and  the  thermomechani¬ 
cal  deformation  problems  are  given.  In  each  case,  a  finite  element  (FE)  discretiza¬ 
tion  of  the  underlying  partial  differential  equations  is  used  and  will  not  be  de¬ 
scribed  in  detail.  Overall,  the  optimization  problem  will  be  to  place  material  in  a 
design  domain  denoted  Q,  which  is  divided  into  a  union  of  N  cells  or  elements  u f 
(U^a f  =  Q),  which  are  quadrilaterals  in  2-D  or  hexahedron  in  3-D.  Bilinear  and 
trilinear  basis  and  testing  functions  are  used  in  2-D  and  3-D,  respectively. 

2.1  Heaviside  Projection 

The  Heaviside  projection  method  (HPM)  for  topology  optimization  was  first  intro¬ 
duced  by  Guest  et  al.1  and  provides  explicit  control  over  minimum  length  scale  in 
addition  to  eliminating  checkerboard  solutions.4  In  the  HPM  (in  contrast  to  sensitiv¬ 
ity  filtering  or  density  filtering)  the  design  space  (</>)  is  separated  from  the  physical 
space  (p)  or  mesh  elements.  Typically,  the  design  variable  locations  are  coincident 
with  the  mesh  nodes  in  physical  space,  though  this  is  often  out  of  convenience — 
the  design  variables  may  exist  anywhere  in  the  domain.  As  will  be  seen,  minimum 
length  scale  control  is  achieved  by  linking  the  physical  space  to  the  design  space 
through  a  neighborhood  with  a  specified  minimum  radius,  rm jn.  In  this  way,  only 
the  design  variables  0*  within  the  specified  minimum  length  scale  will  contribute  to 
the  density  pe  of  element  e. 

More  specifically,  each  element  (which  corresponds  to  physical  variables)  is  as¬ 
signed  a  local  neighborhood  set  containing  design  variables,  defined  as 

ieNl  if  ||xj  -  xe||  <  rmin,  (1) 

where  x,(  is  the  location  of  the  ith  design  variable  and  xe  is  the  location  of  the 
centroid  of  the  eth  element.  This  neighborhood  is  then  used  to  define  an  intermediate 
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variable  pe,  which  is  a  weighted  average  of  design  variables 


e  _  TsieNt  <MX*  -  X6) 

A4  (  ~e\  ’  (2) 

EjeNi  w(xi  -  x  ) 

where  weighting  function  w  is  defined  as 

w(x*  -  xe)  =  1  -  — — — .  (3) 

Aran 

Finally,  the  unpenalized  physical  variables  p  are  defined  via  p  as 

pe  =  H(pe(4>))  =  1  -  *;w.  (4) 

where  Ff  (•)  is  the  regularized  (continuous  approximation)  Heaviside  function  and  f3 
controls  the  curvature  of  the  regularization  and  may  range  from  5  to  upwards  of  50. 
Note  that  as  f3  approaches  infinity,  the  regularized  (continuous)  Heaviside  function 
approaches  the  true  discontinuous  Heaviside  function.  While  the  design  becomes 
more  crisp  (i.e.,  the  width  of  the  transition  regions  from  ff  =  1  to  ff  =  0  is  reduced) 
as  /3  increases,  the  function  becomes  more  nonlinear,  making  the  optimization  more 
susceptible  to  suboptimal  local  minima. 

A  penalization  scheme  is  necessary  to  encourage  convergence  to  0-1  designs  that 
have  minimal  intermediate  values  of  pe.  For  a  given  elemental  density  ff,  a  penal¬ 
ization  reduces  the  element  stiffness  while  leaving  the  density  used  for  computing 
the  volume  fraction  unchanged.  This  reduction  in  stiffness  then  causes  intermediate 
densities  to  be  inefficient.  The  most  common  penalization  scheme  is  solid  isotropic 
material  with  penalization  (SIMP),5,6  given  by 

=  (//)"(!  Anin)  WW  (5) 

where  p  is  a  penalization  exponent  and  pmin  is  a  minimum  density  used  to  avoid 
setting  the  stiffness  to  zero. 
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2.2  Compliance 


For  the  static  compliance  minimization  problem,  the  goal  is  to  optimize  the  distri¬ 
bution  of  material,  again  denoted  by  p,  across  the  design  domain  so  as  to  minimize 
external  work  subject  to  a  volume  of  material  constraint.  Physical  variable  pe  repre¬ 
sents  the  volume  fraction  of  element  e,  which  modifies  the  material  properties  of  the 
element.  The  resulting  optimization  formulation  takes  on  the  following  well-known 
form: 

min  c(<p)  =  FTd 
4> 

subject  to:  K(</>)d  =  F 

N  (6) 

5>(0K<V' 

e=l 

0  <</>'<  0max  =  1  Vi  efl 

where  F  is  the  vector  of  applied  nodal  loads,  d  is  the  vector  of  nodal  displacements, 
K  is  the  global  stiffness  matrix,  ve  is  the  volume  of  element  e,  and  V  is  the  total 
allowable  material  volume.  For  later  use,  let  us  also  define  the  variable  Vf  as  the 
total  allowable  volume  fraction,  or  the  ratio  of  V  to  the  total  volume  of  the  design 
domain. 

2.3  Approximate  Process  Model _ 

A  high-fidelity  AM  process  model  must  consider  numerous  aspects  that  would  be 
difficult  to  fully  capture  in  an  optimization  algorithm.  For  example,  selective  laser 
melting  involves  a  laser  impinging  on  a  powder  bed  and  melting  the  powder  in 
an  area  determined  by  powder  depth,  absorption,  temperature,  laser  power,  focal 
length,  and  beam  width.  After  melting,  the  material  solidifies  and  heat  is  trans- 
fered  to  the  surrounding  material,  including  both  built  material  and  powder.  The 
heat  is  also  transfered  to  the  surrounding  environment  and  build  plate.  While  the 
consolidated  material  cools,  it  undergoes  thermal  deformation  and  may  warp  the 
surrounding  built  material.  This  thermal  deformation  is  the  most  important  aspect 
to  capture  for  this  study,  and  so  we  attempt  to  approximate  this  process  with  the 
surrogate  models  discussed  below. 

In  both  models  described  next,  the  objective  function  is  defined  as 

dtherm  =  P*  ||lle||,  (7) 

e 
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where  ||ue||  is  the  magnitude  of  the  displacement  at  the  centroid  of  element  e. 

2.3.1  Basic  Thermomechanical  Model 

The  first  approach  is  a  basic  thermomechanical  surrogate  model  in  which  linear, 
isotropic  thermal  expansion  is  used.  The  optimization  formulation  is  then  modified 
to 

min  dtherm(</>) 

subject  to.  K(0)d  Ftherm 

N 

X>(<£  )»e  =  V 

e=l 

where  Ftherm  are  thermal  loads  induced  by  the  thermal  stress  (defined  by  thermal 
expansion  coefficient  a  and  temperature  difference  AT)  and  the  bounds  on  </>  are 
also  constrained  as  before.  Note  that  the  total  volume  constraint  is  now  an  equal¬ 
ity  constraint,  though  this  may  be  relaxed  when  the  2  optimization  problems  are 
combined  into  a  compromise  objective  function.  The  build  plate  is  modeled  as  a 
fixed  displacement  boundary  condition  along  one  edge  of  the  design  domain  and 
is  distinct  from  the  boundary  conditions  applied  in  the  compliance  minimization 
problem  described  above.  This  approach  is  similar  to  that  of  Reference  7,  though 
the  thermomechanical  problem  is  solved  separately,  with  different  boundary  condi¬ 
tions. 

2.3.2  Element  Birth  Model 

The  element  birth  model  involves  the  iterative  activation  of  elements  in  an  FE  mesh 
to  simulate  the  AM-like  process  of  adding  solid  material  to  a  structure.  In  this 
model,  the  physics  involved  remains  purely  thermomechanical,  though  the  iterative 
placement  of  material  changes  the  overall  deformation  of  the  resulting  structure.  As 
such,  no  thermal  diffusion  is  used,  though  it  will  be  added  in  the  future.  Instead,  a 
linear  temperature  history  is  assumed  for  each  element,  starting  at  T0  and  decreas¬ 
ing  by  an  amount  AT  at  each  step.  This  temperature  history  is  an  oversimplification 
of  the  process  and  only  used  here  for  simplicity.  In  general,  the  temperature  profile 
of  an  element  should  be  measured  experimentally  and  will  not  be  a  simple  linear 
function  of  time. 

Figure  1  illustrates  the  element  activation  process  and  per-element  temperature  his¬ 
tories.  First,  all  elements  are  initialized  as  inactive.  An  inactive  element  has  its 
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Young’s  modulus  E  set  to  a  low  value  Tmm  and  its  thermal  expansion  coefficient  a 
set  to  zero.  Next,  the  first  element  (or  set  of  elements)  is  activated  and  its  temper¬ 
ature  reduced  by  AT  and  the  thermal  deformation  of  the  entire  mesh  is  computed. 
Note  that  when  elements  (other  than  the  first)  are  activated,  they  may  exhibit  er¬ 
roneous  strain  since  the  element  will  have  been  distorted  by  the  prior  activated 
elements.  This  strain  is  subtracted  from  the  current  strain  when  used  to  compute 
the  stress,  so  that  the  initial  strain  of  the  activated  element  is  zero.  Then,  the  next 
element  is  activated  and  the  process  is  repeated  until  all  elements  are  activated.  For 
this  model,  the  thermomechanical  deformation  (defined  in  Eq.  7)  at  the  last  time 
step  is  used  to  compute  the  objective  function. 


Inactive 


Active 


Deformed  Next  element 

mesh  activated 


Inactive 

6B-H 

T0  /  T0  -  AT 

H-n 

T0  -  2AT  T0  -  3A T 


(a) 


(b) 


Fig.  1  Schematic  of  a)  the  element  activation  process  and  b)  the  element  temperature  histories 


Note  that  the  strain  subtraction  procedure  is  an  approximation,  and  a  more  accurate 
method  would  be  to  reset  the  initial  node  positions  of  the  element,  recompute  the 
element  stiffness  matrix,  and  reassemble  the  global  stiffness  matrix.  It  can  be  shown 
that  subtracting  the  strain  amounts  to  a  zeroth  order  approximation  of  the  updated 
element  stiffness  matrix;  thus,  for  small  displacements  this  approximation  should  be 
sufficient,  though  error  may  accumulate  during  the  build  process  as  the  upper-most 
inactive  elements  deform.  However,  stress  equlibrium  and  energy  conservation  are 
preserved  between  activation  steps. 
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2.4  Compromise  Objective  Function 


The  objective  function  is  then  a  linear  combination  of  the  structural  problem  and 
the  thermomechanical  deformation  problem,  and  can  be  defined  as 

min  /(</>)  =  — c{<t> )  +  - — —  dthenn(<£),  (9) 

<t>  Wc  wd 

where  0  <  w  <  1  is  the  compromise  weight  and  wc  and  wd  are  normalization 
factors.  Finally,  the  widely  used  method  of  moving  asymptotes  (MMA)  is  used  as 
the  optimization  algorithm.8 

3.  Results 

Several  results  of  both  optimization  processes  are  given  in  this  section.  First,  a  set 
of  cantilever  designs  in  a  square  design  domain  are  given  while  moving  the  build 
plate  orientation  to  each  edge  in  the  design  domain.  These  results  are  compared  with 
solutions  generated  using  an  overhand  constraint  approach.  In  addition,  an  example 
using  a  larger,  rectangular  design  domain  is  given. 

3.1  Square  Design  Domain 

In  this  section,  a  square,  lxl  region  is  used  as  the  design  domain  and  4  boundary 
condition  configurations  are  explored,  as  shown  in  Fig.  2.  In  Fig.  2,  orange  dashed 
regions  indicate  the  fixed  boundary  used  in  the  thermomechanical  distortion  prob¬ 
lem  (i.e.,  the  build  plate),  the  blue  dashed  region  indicates  the  fixed  boundary  for 
the  compliance  problem,  and  the  orange  arrow  indicates  the  load.  For  the  pure  com¬ 
pliance  problem,  the  boundary  conditions  result  in  the  cantilever  solution  shown  in 
Fig.  3.  The  4  cases  then  represent  placing  the  thermomechanical  boundary  condi¬ 
tion  on  each  edge  of  the  design  domain. 

The  parameters  used  in  this  problem  were  volume  fraction  Vf  =  0.5,  Heaviside 
regularization  f3  =  10,  penalization  exponent  p  =  4,  compromise  weight  w  =  0.5, 
discretization  size  0.02  (50  elements  per  linear  dimension),  and  filter  radius  rmin  = 
0.1.  The  optimizer,  MMA,  was  run  for  150  iterations.  The  material  properties  for 
all  examples  were  Young’s  modulus  E  =  1  GPa,  Poisson’s  ratio  v  =  0.25,  and 
thermal  expansion  coefficient  a  =  10-5.  Further,  the  temperature  difference  for  the 
thermomechanical  problem  was  AT  =  —  50  K.  Normalization  factors  wc  and  wd 
were  set  to  the  values  of  the  individual  objective  functions  c  and  dtherm  at  the  initial 
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guess.  The  optimization  results  are  shown  in  Fig.  4,  and  the  corresponding  objective 
function  values  are  listed  in  Table  1.  The  objective  function  values  in  Table  1  show 
that  2  results  are  optimal  trade-offs  (designs  b  and  c)  and  2  are  dominated  (a  and  d). 
In  other  words,  the  compliance  and  thermal  displacement  objective  function  values 
of  designs  a  and  d  are  both  strictly  greater  than  those  of  both  designs  b  and  c.  For 
designs  b  and  c,  design  b  has  lower  compliance,  but  higher  thermal  displacement 
than  design  c;  thus,  these  designs  are  optimal  trade-offs  in  that  neither  outperforms 
the  other  in  terms  of  both  objectives. 


(c) 


(d) 


Fig.  2  Schematics  of  the  4  design  cases  with  a  square  design  domain 
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Fig.  3  Static  loading  solution  (i.e.  ,W  =  1) 


Table  1  Objective  function  values  for  the  four  lxl  design  cases 


Problem 

Case  a 

Case  b 

Case  c 

Case  d 

c  (kNm) 

34.9 

25.4 

26.6 

27.3 

dtherm  (l*m) 

825 

390 

236 

519 
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Fig.  4  Results  of  the  4  design  cases 


The  thermomechanical  method  discussed  here  has  a  similar  goal  to  and  may  be 
compared  with  an  overhang  constraint  approach.3  The  overhang  constraint  method 
was  applied  to  the  design  cases  illustrated  in  Fig.  2  using  an  overhang  angle  of 
45°,  and  the  results  are  shown  in  Fig.  5.  (Note  that  Fig.  5  also  shows  a  density  plot 
though  a  different  color  scheme  is  used.)  Two  design  cases  (b  and  c)  appear  fairly 
similar  between  the  2  methods,  while  cases  a  and  d  differ  more  significantly.  The 
advantage  of  the  thermomechanical  approach  is  that  the  overhang  constraint  may  be 
overly  restrictive  for  short  overhanging  sections  of  a  part.  In  other  words,  the  length 
of  an  overhanging  section  is  also  important  in  addition  to  overhang  angle  to  build 
success.  The  overhang  constraint  approach  does  not  consider  the  length  of  a  section, 
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while  in  the  thermomechanical  approach  this  may  be  a  natural  outcome.  Designs  b 
and  c  may  then  be  similar  in  both  methods  because  the  overhang  constraint  does  not 
significantly  alter  the  designs  when  compared  with  the  unconstrained  design  (see 
Fig.  3)  and  so  neither  does  the  thermomechanical  method.  Design  cases  a  and  d  do 
differ  significantly  from  the  unconstrained  design,  and  so  because  the  thermome¬ 
chanical  method  allows  for  overhanging  features,  it  does  result  in  designs  different 
from  the  overhang  constrained  designs. 


(c)  (d) 


Fig.  5  Results  of  the  4  design  cases  using  a  45°  overhang  angle  constraint 


Figures  6  and  7  show  the  thermally  induced  displacement  of  design  cases  a  and  d. 
In  each  figure,  subfigure  a)  shows  the  initial  configuration  and  subfigure  b)  shows 
the  final  configuration  after  cooling.  Note  that  the  displacement  is  exaggerated  by 
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a  factor  of  1000,  and  elements  with  a  density  of  less  than  0.25  are  removed  as 
large  deformations  associated  with  these  elements  obscure  the  results. 


(a) 


(b) 


Fig.  7  Thermal  deformation  of  design  case  d:  a)  initial  geometry  and  b)  final  result 
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To  illustrate  the  effect  of  varying  the  compromise  weight  w,  design  case  d  from 
Fig.  4  was  repeated  using  8  values  of  w:  0.1,  0.2,  0.3,  0.4,  0.6,  0.7,  0.8,  and  0.9. 
These  results  are  shown  in  Fig.  8,  where  Fig.  8a  shows  the  result  for  w  =  0.1,  and 
w  increases  to  a  value  of  0.9  for  Fig.  8h.  (Recall  that  the  result  in  Fig.  4d  corresponds 
to  a  value  of  w  =  0.5.)  As  can  be  seen,  for  lower  values  of  w  the  results  display 
fewer  overhanging  features  and  may  represent  designs  that  are  less  prone  to  build 
failures  with  an  additive  process.  In  addition,  Fig.  9  shows  the  individual  (penalized) 
objective  function  values  (c  and  <7lhemi )  with  the  lowest  value  of  c  corresponding  to 
the  highest  value  of  w. 
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(a)  w  0.1  (b)  w  0.2 


(f)  w  =  0.6 


(g)  w  0.7  (h)  w  0.8 

Fig.  8  Results  of  design  case  4  with  varying  w 
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Fig.  9  Individual  objective  function  values  for  varying  w  in  Fig.  4 


Finally,  design  case  b  was  repeated  using  the  element-birth  model.  In  this  example, 
the  discretization  size  was  reduced  to  30  elements  per  linear  dimension,  and  the 
stiffness  penalization  exponent  was  increased  to  pi  =  6.  The  size  of  the  added  ma¬ 
terial  at  each  step  was  0.1-by-0.1,  or  3  elements,  and  the  temperature  was  decreased 
for  each  active  element  by  0.5  K.  Optimization  results  are  shown  for  compromise 
weight  w  —  0.9  and  w  =  0.95  in  Fig.  10  along  with  the  pure  compliance  result 
(w  =  1)  for  comparison. 


(a)  w  =  0.9  (b)  w  =  0.95  (c)  w  =  1 

Fig.  10  Results  of  the  element-birth  model  with  varying  w 
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3.2  6x1  Design  Domain 


In  this  example,  a  6  x  1  design  domain  was  used  with  half-symmetry  boundary 
conditions.  The  bottom  corners  were  fixed,  and  load  was  applied  to  the  top  center. 
Parameters  for  the  optimization  algorithm  and  discretization  were  identical  to  those 
of  the  1  x  1  design  cases.  First,  the  thermomechanical  fixed  boundary  was  placed 
on  the  bottom  edge  of  the  domain,  resulting  in  the  design  shown  in  Fig.  11a.  Also 
shown  in  Fig.  lib  is  the  same  compliance  problem  optimized  with  an  overhang 
angle  constraint  of  45°.  It  is  apparent  from  Fig.  11  that  the  thermomechanical  re¬ 
sult  still  exhibits  overhanging  features,  though  note  that  overhanging  features  may 
not  lead  to  a  failed  build  because  the  length  of  those  features  is  important.  The 
individual  objective  function  values  for  the  thermomechanical  optimization  were 
c  =  4242.17  N  m  and  dtherm  =  59.2  pm. 

Finally,  the  problem  was  repeated  with  the  thermomechanical  fixed  boundary  on 
the  top  edge  of  the  region.  This  result  is  shown  in  Fig.  12  (also  compared  with 
an  overhang  constrained  design)  and  had  individual  objective  function  values  of 
c  =  4712  N  m  and  dtherm  =  76.8  pm.  In  this  case,  building  from  the  top  down  is 
clearly  suboptimal  as  both  the  compliance  and  thermal  displacement  are  higher. 


(b) 


Fig.  11  Results  of  the  6x1  design  domain  using  a)  the  thermomechanical  model  and  b)  a  45° 
overhang  constraint 
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(a) 


(b) 


Fig.  12  Results  of  the  6x1  design  domain  (with  the  build  plate  placed  on  the  top  edge)  using 
a)  the  thermomechanical  model  and  b)  a  45°  overhang  constraint 


3.3  Comparison 

Because  the  element-birth  approach  is  much  more  computationally  expensive  than 
the  basic  thermomechanical  model,  it  is  worthwhile  to  compare  the  results  of  the 
approaches  directly.  If  the  basic  thermomechanical  optimization  approach  yields  re¬ 
sults  that  are  near-optimal  using  the  element-birth  model,  we  may  use  the  cheaper 
optimization  result  as  an  initial  guess  for  the  more  expensive,  hopefully  saving  com¬ 
putation  time. 

To  that  end,  the  basic  thermomechanical  model  was  optimized  using  the  same  dis¬ 
cretization  as  the  element-birth  model,  and  the  results  are  shown  in  Fig.  13.  After 
optimization,  the  thermal  distortion  objective  (dlherm)  of  the  basic  result  was  evalu¬ 
ated  using  the  element-birth  model.  In  this  case,  the  unpenalized  thermal  distortion 
was  1.78  x  10-4  for  the  element-birth  result  and  2.14  x  10-4  for  the  basic  thermo¬ 
mechanical  result,  which  is  about  a  20%  difference. 

3.4  3-D  Results 

Finally,  the  basic  thermomechanical  model  was  used  in  3-D  in  a  1  x  0.5  x  0.5  region 
in  a  cantilever-like  loading  with  a  volume  fraction  of  Vf  =  0.3.  The  fixed  boundary 
was  located  on  the  x  =  0  plane,  and  the  load  was  — ^/-directed  and  located  at  the 
point  (1, 0, 0.5).  Two  build  directions,  x  and  y — meaning  that  the  build  plates  were 
located  on  the  x  =  0  and  y  =  0  planes,  respectively — were  tested  using  a  weight  of 
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w  =  0.5.  The  results  of  the  optimization  are  shown  in  Fig.  14  with  the  pure  com¬ 
pliance  (w  =  0)  result  shown  in  Fig.  14(a)  for  comparison.  The  compliance  of  each 
result  was  35.4,  40.8,  and  43.4  kNm  for  the  pure  compliance,  x  build  direction,  and 
y  build  direction  results,  respectively.  Comparing  the  thermal  displacement  between 
the  2  build  directions  is  not  informative  because  the  surface  area  of  the  build  plates 
is  different;  however,  comparing  the  thermal  displacement  of  the  pure  compliance 
problem  with  each  reveals  that  each  result  is  a  trade-off  with  the  pure  compliance 
result,  as  expected.  Table  2  lists  the  objective  function  values  for  each  result. 


(a)  (b) 

Fig.  13  Comparison  of  a)  basic  and  b)  element-birth  thermomechanical  optimization  prob¬ 
lems 


Table  2  Objective  function  values  for  the  3-D  design  cases 


Problem 

Pure  compliance 

a? -directed  build 

y  -directed  build 

c  (kNm) 

35.4 

40.8 

43.4 

x -directed  c/therm  (pm) 

593 

411 

652 

-directed  dtherm  (pm) 

443 

380 

238 
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Fig.  14  Results  for  the  1  x  0.5  x  0.5  design  domain  with  a)  no  thermomechanical  model  (w  =  1), 
b)  build  direction  in  the  x-direction  and,  c)  build  direction  in  the  ^/-direction 
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4.  Conclusions 


An  optimization  technique  was  presented  that  pairs  a  minimum  compliance  prob¬ 
lem  with  surrogate  models  for  the  AM  process.  The  goal  was  to  minimize  ther¬ 
mal  distortions  induced  during  the  build  process  to  generate  more  readily  manufac¬ 
turable  designs.  Two  process  models  were  used,  the  first  being  a  one-step  thermo¬ 
mechanical  simulation  in  which  the  entire  part  begins  at  the  same  initial  tempera¬ 
ture  and  cools  at  the  same  rate.  The  second  was  an  element-birth  iterative  model,  in 
which  elements  are  activated  in  turn  and  at  a  given  initial  temperature.  At  each  step, 
activated  elements  are  cooled  at  a  linear  rate.  For  the  optimization  formulation,  an 
HPM  and  SIMP  penalization  were  used  as  the  topology  parameterization,  and  the 
MMA  algorithm  was  used  as  the  optimizer.  Finite  differences  were  used  to  compute 
the  gradient  of  the  topology  representation. 

Results  demonstrate  the  viability  of  the  method,  and  several  examples  were  pre¬ 
sented  with  different  design  domains.  Further,  the  results  were  compared  with  over¬ 
hang  angle  constrained  results,  demonstrating  some  similarities  and  differences 
between  the  methods.  Finally,  the  basic  thermomechanical  model  was  compared 
with  the  element-birth  model  demonstrating  that  the  2  methods  do  give  different 
optimization  results,  even  though  they  use  the  same  objective  function.  The  dif¬ 
ference  in  the  thermal  displacement  objective  function  value  for  the  same  design 
case  was  found  to  be  approximately  20%.  While  this  is  a  rather  large  difference,  the 
cheaper  model  may  be  used  for  an  initial  guess  for  the  more  expensive  model,  as  the 
element-birth  model  is  significantly  more  computationally  expensive  for  a  full  op¬ 
timization  run.  Consider,  the  computational  complexity  of  a  naive  implementation 
of  the  basic  model  is  0(N 4)  at  worst  when  using  a  finite  difference-based  gradi¬ 
ent  computation  because  a  matrix  equation  must  be  solved  for  each  element  in  the 
mesh.  For  the  element-birth  model,  the  complexity  grows  to  0(N5)  because  each 
step  of  the  element-birth  process  involves  solving  a  new  matrix  equation,  which 
must  be  done  for  each  element  again  to  compute  the  gradient.  Regardless  of  the 
matrix  equation  solution  technique  (which  may  reduce  the  computational  complex¬ 
ity),  the  element  birth  model  will  always  have  a  computational  complexity  of  one 
power  higher.  Additional  efficiency  may  of  course  be  gained  by  using  an  adjoint 
gradient  approach  for  the  basic  thermomechanical  model,  though  it  remains  to  be 
seen  if  an  adjoint  can  be  formulated  for  the  element-birth  model.  In  the  future,  the 
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designs  generated  with  both  approaches  will  be  tested  using  a  suitable  AM  tech¬ 
nique  to  determine  buildability. 

Using  the  concept  of  pairing  a  process  model  with  the  optimization  problem,  there 
are  several  other  options  that  may  yield  better  results  and  will  be  explored  in  the 
future.  First,  residual  stress  minimization  could  be  used  rather  than  displacement 
minimization.  Second,  build  failures  often  occur  because  the  built  part  warps  to 
such  a  degree  that  it  is  impacted  by  the  recoater  blade  or  roller.  A  constraint  could 
then  be  placed  on  the  maximum  allowable  thermal  deformation  seen  in  any  single 
build  layer. 
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List  of  Symbols,  Abbreviations,  and  Acronyms 


2- D  2-dimensional 

3- D  3 -dimensional 

AM  additive  manufacturing 
FE  finite  element 
HPM  Heaviside  projection  method 
MMA  method  of  moving  asymptotes 
SIMP  solid  isotropic  material  with  penalization 
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